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\ w;''  Qi^z  S ' 

S abstract  ( jJ/WC  '> 

A model  system  of  equations  which  has  been  used  to  describe  the  miscible 
d isplacement  of  one  incompressible  fluid  by  another  in  a porous  medium  is 
the  coupled  quasi  linear  system^fop  c = c(x,t)  and  p = p(x,t)  for 
x e il,  t t (0  ,T]  given  by  \ 

V*[a(x,c){Vp  - y(x,c\vg}]  -V-u  = f ^ (x)  , 

V • 1 b ( x , t ) Vo  ] - u (x,c,  yp)  *Vc  =v,(x)  — - - f^(x,c)  , 

\ 

with  appropriate  initial  and  Neumann  boundary  conditions.  Another  case 
considered  is  when  b = b(x,c,Vp)  above.  ^ Iterative  methods  are  presented 
and  analyzed  which  are  based  on  using  a preconditioned  conjugate  gradient 
iteration  for  approximately  solving  the  systems  of  linear  equations  produced 
at  each  time  step  by  Galerkin  methods  for  time-stepping  the  above  system. 
Optimal  order  convergence  rates  are  obtained  for  the  iterative  methods.  The 
iterative  methods  are  computationally  more  efficient  than  Galerkin  methods 
previously  proposed  to  solve  the  above  system.  The  use  of  different  time 
increments  in  the  t.ime-steppi  nq  procedures  for  the  different  variables  is  also 
presented  and  analyzed.  The  use  of  unequal  time  increments  takes  advantage  of 
different  smoothnesses  in  time  of  the  physical  variables  p and  c and 
greatly  reduces  the  work  done  in  the  computation  of  the  approximate  solution. 
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SU'.NU' U'/VNi  i:  ANP  CXI'l.ANA  I'U'N 


rho  numerical  approximation  of  a ayat  em  of  partial  different  IaI  oqua- 
t tons  usevt  to  mtHle  1 the  miscible  d ispl  acement  of  one  incompressible  fluid 
by  .mot  her  in  i<otous  modi  a is  considered.  Pol  example,  the  model  system 
tiAS  been  used  to  desci  ibe  the  variables  of  pressui  e Alld  the  chanqinq 
concentration  of  a ehemieAl  solvent  in  o 1 1 used  to  flood  oil  wells  to  push 
the  i>il  through  the  porous  media  toward  product  ion  wells  to  act  ureal  or 
recovery  of  oil  fi.im  unde rq round  reservoirs. 

Wheeler  and  the  author  have  recently  presented  some  numerical  methods 
tor  approximatin'!  t h«*  solution  of  the  model  system  and  have  obtained 
optimal  order  etiot  estimates  for  these  methods.  This  paper  presents  and 
analyzes  methods  which  air  . -output  at  tonal  ly  mote  efficient  than  eailiei 
methods' . 

In  t inie-st  t'ppinii  problems  with  time-dependent  coet  f ic  tent  s , numetic.il 
methods  produce  different  systems  ot  1 meat  e>|uat  ions  which  must  be  solved 
at  each  time  step.  Kactotiii'i  a ditfetent  mattix  at  each  time  step  in  t In- 
solation process  can  be  vet  y .•omput  at.  luna  1 ly  expensive  sine*'  the  matt  ices 
are  frequently  ot  the  order  I, t'lH>  ' l,PPP.  The  methods  presented 
in  this  p.uei  require  the  t act  or  i/at  ion  ot  only  two  matri.es.  Then 


an  iterative  procedure  which  compensates  tor  the  tact  that  new  matrices 
are  not  factored  at  each  time  step  is  utiixi  to  obtain  approximate  soluti. 
Sv. mil  leant  amount  s of  ..'input  at  ion  are  saved  by  these  methods.  optimal 

older  estimates  on  rates  ot  eonvei  >|enc«>  tor  t hose  methods  are  obtained. 

__'v 

I 1 

| Dint 

Ml 


rhe  lesponsiblllty  t >t  the  wo  I .1  1 II<1  and  views  expl  essed  III  tills  desci  ip 
t ivi  summat  l 1 1 " with  Ml..',  and  not  with  t h.  ant  hoi  ot  this  report. 
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I'ROULKMS  IN  POROUS  MKIHA 
Richard  K.  l-viiui 


1.  Introduction 

In  1 10|  numerical  approximations  by  Gaierkin  methods  to  it  probl tin  aiism,)  in  tin* 
miscible  displacement  ot  one  incompressible  fluid  by  another  in  a porous  medium  are 
piesented  and  analysed.  A brief  discussion  of  the  physical  problem  is  iiuvti  in  lid] 
and  a mathemat  ica  l model  which  is  sufficiently  qeneral  to  incot  (km  at  e most  ot  the 
ilia  jot  features  ot  the  physical  problem  (see  [10,l‘>,17))  is  presented.  In  this  paper, 
we  -hall  present  and  analyze  some  methods  for  t ime-st  eppinu  these  nuxiel  equat  ions 
which  aii-  much  more  efficient  comput at ional ly  than  the  methods  discussed  m 110] . We 
•hall  use  an  iterative  method  based  on  a precondi  tioned  conjugate  qiadient  it  elation 
to  approximate  the  solut  ion  ot  the  systems  of  linear  equations  which  at  ise  ttom  a 
Galeikin  approx imat ion  of  the  model  equations.  The  iterative  methods  presented  pre- 
eive  the  accttiaey  mil. -lent  in  the  underlying  Gaierkin  metliod  and  let  us  obtain  vei  y 
neai ly  optimal  possible  orders  for  the  wot k involved  in  solving  the  lineal  systems  ot 
equat  ions . 

We  tust  present  tin  model  equations  for  our  physical  problem.  Kind  c ctx,t) 
and  p p(x,tl,  solutions  of 


(1  . n 

l l • ' 


V*[a(Vp  - qV.il]  -V • u f 
V-  l bVc  I - a -Vo  - * 5‘*  - f , . 


t O I X t , 


«».Tl  . 


Whet  e 


t.'  is  a bounded  domain  in 


with  boundary 


r and  uix.c.Vp)  -alVp  - > V»i ) is  a veetoi  in  K* . Here  a aix,.'.  > i(x,c), 

.1  q(x),  f t ^ (x)  , ^ -Mx!,  t,  - t,(x,c),  and  b are  specified.  Kot  ease  of  expos 

t ion,  we  shall  oonstdei  two  oases  for  b:  I'.iiie  1,  with  b blx,t),  and  Case  II,  with 


b blx.c.p  ,p  1 . 

X y 

i • ib  ise  ot  Case 


The  cases  ot  .neatest  physical  interest  (see  l'.,  1 'I  ate  Case  II  and 

l with  b b(xl  (special  tesults  foi  this  sub,  ase  aie  obtained  in 


Coro  l I ui  y -l . .’  1 . We  i ••  -a  mi.  • that  the  t o I I ow  i nq  hound  at  \ and  mil  i a ! cot  id  i t ions  hold 

• pole-,  »i  e.t  by  III.  United  state  Army  un.let  Coni  i a,  I Numbets  P.VVG.si  '■  C I'i'.'l  m. I 
I 'A  AG.’"  M -C.  i'll. I . This  mateiial  is  based  upon  wot  k supi’oite.l  bv  the  National  s,  i.nc. 
l oundat  n»n  undei  til  ant  Nuinhei  MCs  ,’H-iMSjs , 


i 


(1.3) 


x e 3fi,  t(  .1 


a 


0 , 


(1.4)  b — = 0 , x e 3ft,  t € J , 

a V 

(1.5)  c (x,0)  = cQ(x)  , x e ft  , 

3F 

where  — is  the  normal  derivative  of  F on  the  boundary  of  ft.  Note  that  (1.1)  — (1.5) 

OV 

will  define  p(x,t)  only  to  within  an  arbitrary  constant.  We  shall  normalize  p by 
the  condition  that 


(1.6)  t— j-  / p(x,t)dx  = 1,  t e J , 

m n 

where  |ft|  is  the  measure  of  the  domain  ft. 

We  note  that  the  analysis  that  follows  would  easily  treat  forcing  functions  f 
and  f^  which  are  smoothly  distributed  over  ft.  If  singular  functions  are  used  to 
model  the  effect  of  small  injection  and  production  wells,  the  analysis  will  fail.  Thus 
we  shall  make  the  assumption  that  for  our  problem,  the  sources  and  sinks  are  smoothly 
distributed  and  shall  then,  without  loss  of  generality,  assume  that  f H 0 and 
f^  ^ 0 for  the  remainder  of  the  paper. 

Continuous  time  Galerkin  approximations  for  (1.1) -(1.5)  are  presented  and  analyzed 
in  [10] . By  lagging  or  extrapolating  the  coefficients  in  discrete-time  versions  of 
these  methods,  we  are  able  to  linearize  and  uncouple  the  systems  of  equations  required 
in  the  approximations.  However,  the  discrete-time  versions  of  these  methods  require  that 
different  systems  of  linear  equations  be  solved  at  each  time  step;  this  is  a computa- 
tionally expensive  process.  In  this  paper,  we  shall  present  and  analyze  iterative 
methods  which  require  the  factorization  of  only  two  matrices  for  the  total  solution 
process.  The  use  of  iterative  methods  to  approximate  the  solution  of  the  linear  equations 
arising  from  the  parabolic  equation  is  an  extension  of  the  techniques  developed  in  [6,91 
for  quasi linear  time-dependent  problems.  The  use  of  iterative  methods  to  approximate 
the  solution  of  an  elliptic,  basically  time- independent  equation,  only  to  within  the 
accuracy  of  the  time  truncation  error  from  an  associated  time-dependent  problem  seems 


to  do  new.  We  emphasize  that  unlike  standard  iterative  procedures  for  elliptic  equa- 
tions, all  the  methods  presented  will  only  require  a number  of  iterations  which  is 
independent  of  At,  h,  and  n to  obtain  a norm  reduction  sufficient  to  feed  the 
associated  parabolic  problem.  Although  a preconditioned  conjugate  gradient  iterative 
method  will  be  presented,  any  method  which  achieves  the  specified  norm  reductions  will 
• suffice  in  all  of  the  analysis  to  follow. 

In  the  physical  problem  which  motivates  our  consideration  of  (1.1)  — (1.5) , the 

pressure  p is  much  smoother  in  time  than  the  concentration  c.  In  order  to  take 

advantage  of  this  difference  in  smoothness  of  p and  c,  we  shall  use  different  time 

increments  At  and  At  in  our  analysis  for  the  time-stepping  of  the  systems  of 
c p 

equations  arising  from  the  equations  for  concentration  and  pressure.  By  using  these 
unequal  time  increments,  we  shall  need  to  update  the  pressure  variable  much  less  fre- 
quently than  the  concentration  variable  and  thus  avoid  considerable  unnecessary 
computation . 

In  Section  2 we  introduce  two  families  of  finite  element  spaces  which  we  use  to 
, approximate  our  unknown  functions  p and  c.  We  present  the  hypotheses  on  (1.1)- (1.5) 
and  the  solution  (c,p) , discuss  elliptic  projections  for  c and  for  p,  and  present 
our  basic  Galerkin  approximation  of  (1.1)-(1.5)  together  with  several  modifications. 

Ir.  Section  3 we  present  our  preconditioned  conjugate  gradient  modifications  of  the 
methods  described  in  Section  2 and  analyze  the  effect  of  the  iterative  approximation  on 
a single  time  step.  In  Section  4 we  obtain  global  error  estimates  for  the  various 
methods  described  in  Sections  2 and  3.  Section  5 contains  a brief  discussion  of  the 
estimates  of  the  work  of  computation  for  the  methods  presented  in  this  paper. 
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2 . I t t>Um  1 na r iqs  and  Description  of  Basic  dalerkin  Appro ximat j o its 


Let  (u,v)  - / uvdx,  <u,v>  - / uvdo,  ||u| 
, il  Dil 

K 

W (U)  bo  the  Sobolev  space  on  with  norm 


(u,u) , and  u I * < u#u ).  I < i 


l 2Jfc 


let  Kk  " 3x 


LS(iJ) 


with  the  usual  modification  for  s - <•'.  When  s “ 2,  denote 


VF  - l F , F ^ ) , write  ||Vf||  ^ in  place 


will  denote  the  Sobolev  space  on  3i». 


|F.  ||\  * ||  F J|S  ,1s.  Also  .1*011) 

wK  * w 

s s 


Let  i.V^  •'  family  of  finite-dimensional  subspaces  of  with  the  follow- 


ing property 


For  p “ 2 or  •»,  there  exist  an  integer  r ' 2 and  a constant  K 

0 

such  that,  for  1 «■  q «■  r and  tf>  t Wq(il)  , 

P 

12. .M  inf  (|U  “ \ll  0 ♦ h|W’  “ X II  1 * Ill’ll  h‘l  • 

\«.V'  w w 0 wq 

h P p p 

We  also  define  a family  of  finite-dimensional  subS}vaces  of  H*  (il)  called  (.V^  1 which 

satisfies  the  same  property  as  { } with  r replaced  by  s.  We  also  assume  that 

the  families  (.V^l  and  l.V  1 satisfy  the  following  so-called  "inverse  hypotheses": 

it  » f .V,  and  4’  f ,V.  , 
h h 


L (il) 


" *0h  1 1 W 


b)  ||  V*  || 


- K h 

•v  — 0 

i (in 


IN' II J 1 Koh  !|v*  ||  . 

Best t ict  i!  as  follows  (with  (SI  denoting  the  collection  of  restrict  ions!  : 


iS! : 1.  il  is  M* -regular,  i.e.,  if 


-Av  * ttv  - C , x » 1.  , i Ol  ) , 

Hv 

- n . x r at)  . 

and  ('. , 1 ) *■  < n,  1 ) * 0 it  6 = 0, 

,ht'n  IMI,  • K(»:)(|k||  ♦ II nil  }» 

H Oil) 

• *>:  is  l.ipschitz. 

IV>i  the  follow!  1)4  assumptions,  we  shall  restrict  tin-  variable  >|,  to  u-  between 
■'  physically  determined  constants,  i.e. 

•-1'  "M,  <_  - M , M # > 0 , M*  - 1 , 

* and  M 1 ♦ t , tor  some  t s 0 (set'  (101))  • Assume  the  following 

in  1 ar it  y tot  a,  >,  b,  u,  and  <f  : 

1.  There  exist  uniform  constants  such  that 

* 

d v a,  ■ a (x , 4 j ) - a v Kj  • 

b ' 0 - # »'  (xl  - K ^ , 

')  IXK-Vl  < K)  . 

d>  M 1 Kj  , 


e) 

0 v b v b (x , l 

: ' (v 

•ase  I ) 

t 

f ) 

I'1,  Uol^'l  ,)  i 

' (1  ♦ 

l'b,l>  < 

l 1 

, , q 

> * 

Let 

t he  de 

t IVvit  IVt'ii  ot’ 

a a ( x , c ) 

, b 1 

o(x,t),  > 

> (x.. 

1 , an, 

u 

i 

u (x, 
i 

i ,q)  satisfy 

the  follow i no  as: 

sumpt  i oils  : 

t ot 

i 1 , 

3. 


In  Cast*  II,  when  b b(x,c,Vp)  b (x.c.p^  ) , for  q?,q3  * K< 


(2.7) 


a)  0 " b4  2 , 

b)  b(x,q  ,p  ,p  ) ^ M , 

I x y 

, 3b  it  3b  I t 3b 

O (x,q1,q2,q3)  | + |—  (x , q3  ,q2  ,q  ^ ) | + (x.q^q^q^ 


♦ |~  (x.q  ,q.,,q  ) | * M 
3c  “ ‘ 


Def ine 

(2.8) 


W‘(  (a.b)  ; X ) 
P 


|| 

X Wq(a.b) 
P 


1 < p,  q 


Let  (p,c) , the  solution  of  (1.1)— (l.S),  satisfy  the  following  regularity  assumptions: 

3 c 


(R)  : 


a) 


1.  (JsH  ) 


3t 


i K->  * 


1.  (JjH  ) 


b) 


L (J;H  ) 


3£ 

3t 


K,  * 


I.  (J;H  ) 


(2.9) 


c) 


L (J}H  ) 


a'  2 ♦ t* 

L (J:H  ' 


1 K, 


for  some  f 0 


L (J}W3) 


d) 


L |J|H  ) 


3t 


L ( J;  W 3> 


e) 


32P 


3t 


3 i 

3t 


1 K > 


.1 , 


L (J{H  ) I.  (.1 : H ) 

In  our  analysis,  we  shall  use  a couple  of  auxiliary  elliptic  problems.  This 


technique  was  used  by  Wheeler  in  1 1 1 . Lot  p e .Y  be  the  elliptic  projection  of  p 
into  defined  by 


(2.10) 


( a ( • , c ( • , t ) ) Vp , Yv ) ( a ( • , e ( • , t ) ) Vp , Vv ) 


( a ( • , c ( • , t ) ) ^ ( • , e ( t ) ) Vq , Vv ) , v e A’j  , 


for  each 


t f 0 , where 


U’.U) 


ip(x,t ) 


p(x,t)hix  0.  for  each  t * 


md  where  ({  , c)  is  t lie  solut  ion  of  (1.1)- (1.5).  The  restrictions  (S)  rm;ly  the 
t e 1 low  i no  resul  t [ . 1 8 1 . 

L > n»n.i  .'.1:  1'iier.  exist',  a I'onst  ant  K - K ( ,,i.,K  «K  »K_)  such  that  for  each 
- 3 3 0 1 - 

t * J, 


For  Case  II,  let  \ s 0 be  taken  sufficiently  larue  that  tin  bilineat  toim 


B(*  ,*,) 

(b (c , Vp) V, 

. .Vy.,)  + (u(c. 

Vp)  -Yii  j ,C  ,1  * X (t 

) 

1 

oercive  over  II*  ( ). 

Then  let 

c f be  tiie 

tl 

el  1 ipt ic  ( i o ieot 

ion  ot  c 

lilt 

let  i tied  by 

. . . 1 ))  B ( c , w ' 

B (c,w) 

-(v-  £ . wi  ♦ 

\ ( c , w ) , w « 

* l 

'h  ’ 

■ each  t t . For  Cas 

e I , where 

b l'  (x,  t ) , 

t he  ooef  t'  i c i en  t s 

b (c , Yp' 

i n 

■ . ( liiioii  by  bit)  . Then,  an  in  17,18]  we  can  obtain  the  to  1 low  uni  lemma. 


lemma  f.J:  Then-  exists  a constant  K(  K^li.',b#,\,K  ^ , K ^ , K , 


M)  sncii  t hat  t oi 


or  l , 


a ) c - c i , 

L~  (J  til' ) 


) le  - o)  ! 


. 


I.  ( J ; li  ) 


r It' 

1.  (ft ; H i 

L (vl;H 


b>  !i'-  - f|l  ,v  : ; Klh'  kil  0 


l.  (.) ; L ' 


1,  U1;H  1 


■'  also  make  the  assumpt  ions  on  i.V.  ),  i.Y,  1,  c,  ami  p that  there  exist  a 

h h 

ons*  it  K,  such  t iiat 


l^'i!  , ... 

L (.1 ; t,  ) 


* IN'- II  , , * Ill'll 

iW)  1 CJ,V  ' ,.W) 


Sufficient  conditions  ioi  (2.1‘j)  can  be  found,  as  in  [t>]  and  I’.SJ  . finally,  a m 
I2,i  ,7),  we  can  obtain  the  following  lemma. 

Lemma  2.  l:  Then-  exists  a constant  K , depending  upon  K , K and  K,  such 

t>  0 12 

that 


(2.16) 


3 c 

> 

3t" 


.1, 


3t  ~ ! 


.1 , 


L ( 0 ; H ) L ( J ; H * ) 

The  continuous  time  approximation  of  (1.1)  — (1.5)  car.  be  defined  as  follow:-:  let 
1'  : 10, T)  * .Vh>  the  approx imat  ion  for  p,  and  C : [0,T]  * ,V  , the  approx  imat  ion 
for  c,  be  defined  by  (supressinu  the  dependence  of  the  coefficients  on  x) 


(2.17) 
and 
(2. IS) 
with 
(2. Id) 


(a(C)Vr.Vv)  (a(C)y (C)Vq.Vv) , 


(b(t)VC,Yw)  ♦ (u(C,YP)  *VC,w)  = - (y*  ~ , w)  , 


(*(-)C(-,  0),w)  = (y-  ( • )e(  • .0)  ,w)  , 


w t y , 
h 


w f W , 
h 


where  c is  the  elliptic  projection  of  c defined  in  (2.13).  In  [10]  it  was  shown 
that  the  appropriate  assumptions  from  (K)  and  lyj  yield: 

Case  I : (b  = b [x  , t ) ) 


(2.20)  h 1 1 V ( P p)  ||  ^ , ♦ ||C  - c ||  v ♦ h ||  V (C  - c)||  , 

L ( J ; L~)  L (J;L  ) L^J.-L') 


Case  II:  (b  = b(x,c,Vp)) 
(2.21)  h ' ; V ( P - p)  || 


r s , r ♦ s 
, 1 K?(h  * h ♦ h 


n 


2 + 11°  - C|j  a,  , + h 1 1 V ( C - c)  || 

L (J;L  ) L (J;L“)  L~(J;L~) 


K (h1  + ll^  ' ! 


Let  At  "*  0 , N - T At  f IT,  and  tl  = cult,  o f K . Also  let  ; n(x)  - ix,i 


and  d , 


n t-1  n,  ,, 

(V  - i ) /At . 


In  [10],  the  following  discrete  time  approximation  with  discretization  enor  of 
the  order  At  was  analyzed  and  an  extrapolated  coefficient  Crank-Nieolson-dalerkin 
scheme  was  stated  to  have  \.\t>“  time  d iscret  i zat  ion  error.  Denote  the  aptu  ox imat ion 


of  p by  W : ( 0 t . t 

0 l 


N 


and  the  uperox imat ion  of  c bv 


u-°.\)  to0C -)  .x> . 


II  * II  * 1 ,11*1  II  * 1 

) vw  , V v ! (a(:  )>  u V’ 


.>  i • that  the  I'ot't  f i.  it’nt  in  it 


• i \ .u  i s l n>i  t i out  the  abiebraie  sy:  ■ tem  (.  ..  w;!!i 


Attune t t i e . Ik'Wi'Vi  i , in  m 


.ins’  problems,  t!u  transport  ttim  >;■  laige  oompared 


t.  tin  ilittus'.on  term  and  it  may  be  numerically  advant  .uieous  to  use  with 

vi  v.  it.  thi  iVi't  t i.it'nt  matrix  is  no  1 oner  or  symmetric.  t'hc  lenaindei  .'t  t his 

. i . w i ons  idol  the  ease  j 0. 


!■.  , w"  ‘ does  not  appeal  m (2.22),  we  can  separate  a no.  -A 

. :..t  . Ivina  \ 2 . 22)  at  ttm<  t”*1  and  usuw  that  solution  in  the  oooitieients  tor 
• of  t J . 4 ' at  tin  time  t *'  In  this  way  we  liave  uneoupled.  i ' and 


l ' and  now  must 


on  In  solve  two  separate  linear  systems.  This  sve.it  ly  i edu.es  tin 


,i  ...  ,.i  -.'tot’  lem  and,  con  espond  i no  ly , tin1  work  needed  to  obtain  a soP.it  ion. 

In  tin'  physical  piobleir.  whieh  motivates  out  eons  idol  at  ion  of  ll.l'  il.''<  the 

n*  1 

. • ute  t is  much  smoother  in  t ime  than  t he  eoncent  rat  ion  c . .aus  ''  . i om 

. does  not  differ  radically  from  the  w”  dot ermi rat'd  at  the  previous  tune  level. 
Vi. is  me'ivates  vat  ions  mod  i t i oat  ions  of  l*  1.' ..'4'  whole  equation  i-.'-l  is  solved 

■nl\  at  every  k " t inn1  step  with  k determined  from  the  relative  smoot hnesses 


ime  ot  P 


d . We  shall  eons  niet  several  of  these  modi t icat ions . 


1 ,.f  At,  k.Nt  with  k to  be  ohosen  later  in  different  wavs.  Pons  t del  At  ^ t mu 
levels,  whieh  coincide  with  At -time  levels  but  have  separ.it  ion  At  ( instead  ot  at. 

: ,.t  l*  V u ,p)  be  a part  ieulat  ratio  ot  the  norms  wbioli  enter  into  t ..e  t r unoat  ion 


. .at  tous  del  lvat  ivi's 


detail  • -I  eciii.  samples  latei.' 


and  p.  i We  shall  describe  this  ratio  in  mors- 


(2.25) 
Then 

(2.26) 


One  modification  of  (2 . 22) - ( 2 . 24 ) is  to  choose  k such  that 

k 


* -1 
(0  ) 


At^D  = kAtD  =»  At  . 


In  this  case,  letting  Bn/k  ])  be  the  greatest  integer  less  than  or  equal  to  n/k , we 
replace  (2.22)  (e.g.  in  Case  II  where  b = b(x,c,Vp))  by 

n+1  n 


(2.27)  (y 


- Z 


At 


, x)  * (b(zn,vwIn/kI)vzn+1,vX)  = -(u(zn,vwIn/kl 


) -vz  , X ) , X c «h  . 


Then  W is  determined  by  (2.24)  only  at  the  At^  time  levels.  The  errors  made  in 


the  coefficients  by  lagging  the  pressure  in  this  fashion  is  0(At^D  ) = O(At) . 


Another  modification  which  lias  the  Scune  order  truncation  error,  suggested  by 

3 r„ , [n/k  ] 1 [n/k  1-1 


Todd  Duixsnt,  is  to  replace  VW^I1//k^  by  VKW^n^k 


- VW 


- VW 
2 


i n the 

coefficients  of  (2.27).  This  linear  extrapolation  to  the  midpoint  of  the  current 


At^-time  interval  would  be  as  easy  as  (2.27)  to  implement  (once  two  values  at  At^-time 


levels  have  been  determined)  and  probably  more  accurate  (especially  for  At-time 


levels  near  the  upper  end  of  the  At^-time  interval). 


A more  accurate  modification  of  (2.22)  can  be  defined  by  evaluating  the  coefficients 


b and  u in  (2.27)  at  appropriate  linear  extrapolations  of  the  Vw  from  the  two 


n+1 


previous  At^-time  levels  to  the  At-time  level  t . For  example,  we  can  write 


n = Hn/k]  + v/k  for  some  v = 0,1,..., k - 1.  Then  for  n % k and  v defined  above. 


we  define  to  t>e  the  linear  extrapolation  for  the  time  level  tn  ' from  the 


values  of  F at  the  two  previous  At^  time  leve 


Thus  for  v = 0,1,... ,k  - 1, 


and  0 = (v  + l)/k,  we  define 


(2.28) 


,1  * 8)F,"/k"  - OF1"7"'-' 


r.  ^ k 

n < k . 


Then  i f 


d2F 


< K,  we  are  making  an  error  of  o((At  ) ) by  approximating 


dt 


L ( n , t ; I . ) 


' [ 


- , 0- 


K ( t ‘ ) by  & r . 'Hus  vlrai  ly  defines  k different  ext  ra|H^]at  ions  .r  to  b 

K K 

tepeated  in  ii|ii  iii  c I'riwfcn  t wo  of  the  At  -t  inu*  levels. 

Its  i im  this  ext rapolat  ion,  it  we  choose  V such  that 

_ 1 

* > 

(2.:.'))  k * (At  ft  1 

t hen  wo  see  t hat 

> * i ) * 

to ) ( \t  , )‘  li  k‘  (At  )"U  * At  . 


Cleat  ly  (2.2'1)  allows  k to  he  chosen  quite  1 arqe  and  therefore  W must  be  computed 
much  less  t leniently  than  before. 

With  the  notation  of  (2.28),  oui  new  modification  of  (2.22)  (in  Case  1 when 


b Vi  ( x , t ) i 

can  Vie 

desci  ibed 

as 

(2.11)  (vM(. 

3n.X>  ♦ 

nt  1 

.V'x) 

- (VI  ( 

1 n Case  1 1 , 

1 epl  a.  e 

by 

b ( 2." , 

V n 

; kvw  , 

n .n. 


h 


n ■ k . 


at  the  At  -time  levels  by  (2.24). 

We  also  considet  another  mod  i 1 i cat  ion  of  both  (2.27)  and  (2.11).  We  nott'd  in 

(C)  that  it  the  at'iument  cot  respond inq  to  o in  out  coefficients  satisfies  (2.4),  then 

the  ooeftieieuts  .satisfy  the  bounds  in  (y').  Since  c satisfies  0 c • I (and 

thus  (2.4))  and  line  : 11  is  an  appi  ox  i mat  ion  to  c",  when  evaluatmu  the  coet  t i c i ent  s , 

it  2, " 0,  we  leplace  the  at  qumen  t for  2."  by  0 and  it  2."  s 1 , we  replace  the 

a\  aliment  lot  by  1.  This  type  ot  truncation  of  nrauinents  has  been  discussed 

eatliet  in  14).  rbe  analysis  ot  the  eitoi  that  this  truncation  causes  will  be  contained 

, . * II  , llll  k J 

in  the  proofs  ot  out  mu  joi  results.  We  shall  use  the  notation  u ( , >W  ) m 

(2.27)  to  note  that  the  argument  of  u has  been  truncated  to  lie  between  0 

and  I.  In  the  same  mantlet,  the  modification  ot  (2.11)  with  1 1 uncut. ed  coetticients 


l . W 1 itten  as 


(.  . 


n n • I . . ii  i 1 , . *n  ,ti,  , . n 

(Vd  2.  ,\l  1 (b(t  ' V 2 ,\\i  ~ (a  (2  ,:  VW  ) -V.  ,\1 


\ 


Ive  ii1  ’ t e that  i'i  ,in(  Nio'l  ;oii-sah’l  kill  tin  • t hod-  can  be  del  i tied  as  a not  lie  I iu< 'd  i t i ea  t ion 
o I (2.12!  which  will  have  O ( ( \t  1 ' 1 t I im  t i nival  mil  er  roi  s . See  1 e , " , 1 0 , le  ) tot 
111*  t i.  >ds  ■ > ! till'  li'ttli. 

- I I - 


'•  1 mat  e Solution  ot  t he  Luicat  Kquut  ions  by  Itetation 

tn  thi:.  soot  ion  w>  shall  present  th«>  line. it  equal  ions  art  sum  1 1 .ini  i’. 


tho  i ot  i espotui  i nq  l.’.J4l.  (Wo  note  that  hy  t <»  j « l ao  i nq  f'  VW*1  hy  yw 


(In  k| 


coot  1 ioi.uit  s oi  b and  u in  what  fol  low::  wi  1 l yield  the  lineat  equations 
t rom  t.'.'M.t  We  also  present  an  iterative  method  tot  approximat  i tut  then 

Til-  ooniuqate  qradient  proeedure  presented  here  provides  only  one  ex.unj 
th,  possible  mod.  Moat  ions  ot  U.  O)  and  (.'..'41  that  tall  under  the  analysis 
th-  next  seot  ion.  Any  method  which  provides  the  norm  reduction  defined  in  t 


"ill  preserve  the  results  of  beet  ion  4. 

M. 


1 he  a basis  tot  .V  and  it/..)  2 


i 1 


he  a basis  tot  A'  . We 


l 


denote  the  solution  of  1.')  and  (J.,'41  by  (y. . ,w‘”)  , where 

M. 

a)  s'" 


1 


t t.  1 1 


r ill 

} . U 

» 1 ‘ 1 


M 


b)  W 


y 

..  10 it  > lot  in  such  that  tu/k.  e 

t 1 


We  next  define  several  matt  ices  and  vectors.  Pot  base  11  let 


•l)  ^ UV  . ) ( ivn  . ,u  ) ) , 

M II 


b)  H*11  (0,o)  (bm 


((»((  V f n'\y  ,V„ 

i’  1 k « 1 1 V 


, V)|  > ) 

\ l 


. m . . . m 

r)  V (0  , ) (u  (0 , iO  ) 
\ \ 


( (u  U 


V „"V  i V 


V I 


0 ,11 
l l 


»'  .Vi  ' ■ V f i . , ii  . ) ) 

ll  II 


l 


.HI  . m i.  m * 

A '0'  (■'  (I'll  ( (a  ( ( ' I1  ll  , 1 I V|j  , \‘e  , , 

''T  ' ' i'  a 


t'  1 


, ni  m 

( i'  1 tv  (O') 


* e in  ■* 

' v ( ( ' 0 .it  ,)  l\o,\v 


l I 


t ■)  and 
in  the 

at  isi  n. | 
;ol  ut  ion 
■1  e ot 

0 1 Veil 

ll  I see 

t hen 


e) 


lta(( 


M 

1 

V 


for  i,j  - 1 , . . . , M and  a,  6 = 1,...,M,.  The  matrices  for  Cast:  I ate  correspondingly 

e chosen  in  a very  arbitrary  way.  A Mood  choice 


■impler.  Mete  b and  a can  be 

0 n 


0 


m itjht  be  b - ...  - 

0 0 0 


b(x,c  (x) , VW*  ) and  a 
Vp  ate  mote  or  less  known,  use  these  values  in 
We  can  wtite  (2.32)  and  (2.24)  in  the  form 


a(x,c^(x))  or,  if  average  values  c and 


the  coefficients  fo  leplace 


and  VW 


0 


(3.3) 


i"(U((.n+1  - r.n)  (*  * AtBn (£ ,u>) ) (f.n+1  - c") 
-AtB1'  (f, ,:,')  [,n  + AtO1'  (£  ,(d)  f ",  n / 1 


.m  ...  m _m 
A (£)(u  = T (£) 


and  for  m such  that  m/k  f 22 , 

(3.4) 

We  shall  not  solve  (3.3)  and  (3.4)  exactly;  instead  we  shall  use  a pro-determined 
number  of  preconditioned  conjugate  gradient  I 1 H ,') ! iterations  to  advance  the 

solution  one  time  step.  The  preconditioning  matrices  will  be  chosen  to  be  independent 
of  n.  Specifically,  we  shall  use 


( <.'') 


1,  4>  + AtB 

0 0 


and  lot  the  preconditioners  for  (3.3)  and  (3.4)  respectively. 


Denote  by 


( i.e) 


a)  C 


b)  P 


y amn . and 
i=l  1 1 


y tfnV . , where  m/k  f 5Z  , 
i=l  1 * 


the  approx imations  to  2 and  w'",  respectively,  produced  by  only  approximately  solving 
(t.t)  and  (3.4).  A starting  procedure  for  obtaining  and  will  be  discussed 

later.  Assuming  that  these  quantities  are  known,  we  shall  find  a 


nt  I 


n*  1 


) using  a preconditioned  conjugate  gradient  iteration  to  npproximntt 


(and  thus 

n+1  . n 


. nt  1 


n 


from  ( 1.3) . We  shall  use  different  initial  guesses  for  t."  1 - r 1 1 tot  n 0 and 
lor  n 1.  Wo  shall  use  linear  extrapolation  for  n s 1.  Spoei  f i cal  ly , we  initiali.c 
our  it  elation  for  (3.3)  as  follow?;: 


-1  l- 


At U (a,0)a 


I . (u , 0) x t At  1*  ( u , t' ) a 
0 


Then,  usinq  tlu>  i ni  t ia  l izut  ion  x , q , and  s from  (1.7),  for  i 1,2,...,. 

0 0 0 | 


whore  the  number  of  iterations  ►.  ( will  be  chosen  later,  independently  of  n,  set 


a)  x 


X t 0 S . , whet  I 
1 1 J ) 


- (I.  q , ,q  . ) 

. o i it* 

' " («  ,Ln(«,8UH) 

I I e 


(3.8) 


Id  q 


j ♦ 1 


q . ♦ 0 | . 1.  (a,  8)  s . , 


-1 


- 1 


c)  s.  = h q,  * 0,  . s . , where  0 

jU  0 ‘j  + 1 2 1 i 2 j 


(Lo  qj i l , lj i 1 

%S'V. 


where  (•,•)  is  the  Kuolidean  inner  product.  finally,  set 

ml 


n 

a ♦ x 


l 


In  a similar  fashion,  we  define  the  iterative  apfiroximat ion  ot  (1.41  as  follow: 


a) 


, - k .0 

m - 1:  x0  X0  t>  , 


-2k 

X X 

0 0 


2*k 


(3.10) 


b) 


m 2 0:  q 


q, 


(m* l ) k 


A s A s 
0 0 0 0 


(m  i l 1 k 


mk  ..mk 

A (a)x  - 1 (a)  . 

0 


Usinq  (1.10)  as  an  initial  ization  fen  x , q , and  s , tot  i 1,2,...,*.  - I, 

0 0 0 2 


where  the  numbet  of  iterat  ions  i i will  be  chosen  later,  independently  oi  iti  and 


do tel  mi ne  x . , q . . , and 

] ♦ i V ' 


from  (1.8)  with  1.  and  I.  (it)  evei  ywlien 
1*1  0 


mk 


replaced  by  A and  A (a)  respectively.  Finally  set 

(1.11) 


tm*  Ilk 

il  x 


Clearly  sinci  the  two  iterations  are  intei related. the  order 


ot  computation  must  oo 


1 

• , 1 

2 3 

i . oi  , . 

k 

. . # u t 

k 

* 

k ♦ 1 2k 

> a , ft 

2k  2k+l 

..... 

We  define 

a 

1 

to  be  the  solution 

of  (3. 

» n 

3 ) w i th  f. 

repl  ai 

1 1 n 
i d by  a , 

1 . v . 

-n*  1 

let 

a 

sat  t st  y 

( «.l 

2) 

. n , . ,-n<  l ti 

1.  (a)  (a  - ci 

) -At  B 

n(a,0) 

n 

i * ; 

dp"(a. 

..  n 
n ) a , 

n s |i  . 

o uni 

larly. 

We  do 

1 i m 

,mk 

’ P to  be 

the  so  1 ut  i on  ot 

( i.  1 

3) 

mk 

A 

(a)Pmk 

mk 

1 (a) 

, 

III  > 

1 . 

It  l! 

S Well 

known 

l 1 , 

".Hi  that  thei 

It  exist 

const 

ant  s 

0 < 

Pi  s 1 

and  0 - 

P2  < 

such 

t hat 

a) 

I 

. 0 

1 l 

1 

al* Ci‘  - || 

e Pl  II 

1 

0 2 
I.  (at) 

(it1  - 

o II 

a >11,. 

/ 

b) 

l 

1 

(a) 

<SnM-«nM)|| 

e ■ -',H 

1 

LU  (a)  •' 

(anH 

t a-1' 

'll  , 

n 

n . i ■; 

1) 

l 

| 

l' 

c) 

II 

a\ 

a)J(Pk  - t;k)  ||( 

o ■ ^11 

I 

Ak  (a  1 

(Bk  - 

e°)  || 

e 

' 

n\ 

d) 

l|A,,,k 

1 

(a)*’ 

(H  - 0 1 || 

. 

1 

mk  2 

A ( 1 1 

<Pmk 

, ,(rl)k 
- 2 !*•  ♦ 

k ip, 

m 

whei  e 

the  S 

ubsot  i 

I’t 

e i nd ica t es 

the  Km  1 

idean  norm  of  t ho 

vect oi . 

Given  a 

0 

and 

till  Ic 

exist 

const 

ant  • 

‘ v *,■  *„• 

and  i| 

1 such  that 

M. 


( t. IS) 


T,  tt 

, x L ( t)  x i 

•'>  0 • *0  ~~  1 * . 0 * x t K 1 . 

X I.  X 
1.1 


T mk 


b)  0 v * < V A --&>*  • 


0 T 

y A0y 


If1,  . o * y C IK  ' , 


where  the  constants  are  independent  o.  h and  depend  only  on  the  bounds  tor  the 
ooef  f leients  in  r.ettinq 


1 


a)  <?, 


1 l.  In) 


1 - %/V 


1 + > 


2 

T ' 
2 


b)  Q„ 


1 " 


1 + ( U'q/U'  j ) ■ 


1 ' 

2 


we  know  f iom  ll,‘>,n,8]  that  p . v 2Q . , i = 1 , 2 . It  6 > 0 and 


( l. 17) 


*i  - 6 loq  At /loq  ^7  ' 


i = 1 ,2  , 


then 
( < . I H ) 


P.  < 2 (At)  , 

l 


i = 1,2 


No  1 1 ■ that 


M 


a)  5nM  = }'  an  p. 

i-1  1 1 


(.1.19)  and 


M . 


. . -mk  r -mk 

b)  P = l 6 . . 

i=  1 1 1 


satisfy,  respectively,  for  Case  1, 


a)  („  s".»)  * »,(tn41)v:-"M.vx.  -(a(C*n.V.;'r"l.Vc-".x. 


\ 1 


( 3 . 20 ) and,  for  (n*l)/k  f 72.  , 


b)  (a(C*nU)VfnU,Vy)  = (a(C*nM)>  (C*n* 1 ) Vq.Vy)  , 


V f 


We  note  that  for  n = 0,1,..., N we  have 


(1.21) 


C "(x) 


0 , if  C(x)  v o , 


Cn (x)  , if  0 < Cn (x)  - 1 , 


1 , if  C (x)  ' 1 , 


-■“-a*-1 — 


-In- 


so  that  , suvo  0 ■ 


• * 1, 


n *11, 
- k 


i n n, 
k-  - C 


1 

At  this  j\>  i lit  wr  shall  victim’  some  special  norms  ami  semi -net  ms  t.  om 

let 


mi:  tv- \ . x > . 


b)  I!  x ii  * 


m l 


(hit  )V\,Vy),  for  ease*  1 , 


n | n v n 

((b(o  ,.v  vr  )Vx.Vx). 


for  ease  11  , 


•'  X (a  (t'*')  V\  , V\)  . 


:i  . 


and  ||  • ||  and  1 - • ||  ,u  « un  i • 
a b 


by  (.'.<),  ’ll,  t s «'du  lva  1 ont  to 

V 

, .(in  valent  to  !|v.  !|  toi  a l 1 n. 

We  note  that  in  terms  of  this  not  m-not  at  ion,  (d . 141  and  the  triamiW  uu 


V l e Id 


a'  d - 0*  |,  . ♦ (At  ' ’ ||i'1  - o' 


; i!i«v°  . (a,  )2||«5c°ii  0i  , 


h'  - on“!L  1 (AtM\-,,M  - pn*1|| 


• i l||«Vn|i.  « (At)-  II  A* 


l!rk  - rk 


k ■ - '"I!  k - 

a a 


n ;rmk  . ,’mk, 


I,  rmk  (m-  1 )k  (m-.''  k 


ana  1 ys  i 


>i  n : v 


•qua  lit' 


1 , 


l 


1'ho  ivnvi'r'icni'f  results  ot  Section  4 depend  only  u(*m  the  norm  reductions  (3.24) 


and  not  ujv>n  the  particular  iterative  method  used  to  achieve  these  norm  reductions. 

We  shall  now  define  .1  st. iif  imi  procedure  to  obtain  and  r^.  First  we  must 

0 

compute  ,i  0 which  satisfies 


i »..M 


j - K((h  + (At)  * ) 


One  way  to  obtain  such  an  estimate  is  to  factor  one  additional  matrix  and  solve  (2.13) 

0 . 

dnectly.  Once  P is  obtained,  we  can  use  the  factored  Ai  from  (3.2.g),  a start- 
mu  guess  ot  x(  >'  (or  anything  closer)  in  (3. 10. a)  and  iterate  the  conjugate 
gradient  procedure  «.  times  where  ►.  satisfies  (3.17)  with  3 = 1 to  obtain 


1 ).27) 


Ii'°  - p°I 


Ku1(At)  . 


Komar k : Although  ( < . will  retain  the  proper  convergence  rates  through  all 
the  analysis  to  follow,  foj  bettei  practical  results,  one  should  iterate  k times 
where  «.  satisfies  (3.17)  with  3 2 to  obtain 


l.'°  - P°||  0 , X (At) 
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4.  A ft  tot  1 itiot  list  t mates 

. , ,n  n ,it  n 

lit  this  sivt  ion  we  obtain  t priori  bounds  for  the  erroi  , i - - and  i - p 

tot  the  pi  vs  defined  in  Section  t.  The  first  result.  Theorem  4.1.  states  that 

0 0 1 k , . , 

liven  u startled  procedure  t v> t obtaining  0,1  , 0 , and  r , we  obtain  optimal 

ordet  11  and  l.‘  error  estimates  if  the  iterative  process  reduces  the  error  in 

the  -elution  ot  the  aluebrnie  problems  by  a fixed  factor  (independent  of  n,  h,  and 

At).  Thus  we  need  only  -i  fixed  number  ot  iterations  ot  both  the  concent  rat  ion  and 

pies  -uie  equat ions  at  the  appropriate  time  stops.  We  emphasize  that  the  elliptii 

, th  . , 

j roblom  must  be  approximately  solved  only  -it  every  k time  step  where  k l s 
determined  by  t.'.A'''  and  the  relative  smoothnesses  of  c and  p.  Corollary  4..' 
extends  the  optimal  order  error  results  of  Theorem  4.1  to  the  physioul  - use  whete 
bix'  depends  solely  on  x by  lopluoinq  the  Herat  ion  in  the  equations  foi  the 
'tioontrat  ion  by  a simple  baoksolve  at  each  time  step.  Theorem  -1 . t , the  third  lesult, 
■fates  that,  m Case  II  where  b b(x,e,\p),  by  it  ei.it  inn 

»-  loq  (At  1 /loq  (y'j  1 

t imes  -'tt  t he  equal  ions  for  oonoontiat  ion  at  each  t into  step  and  a t ixed  numbei  ot  t imes 

tot  the  oquut ions  toi  pressure,  we  can  still  obtain  quasi -opt  inn  l I.  error  estimates. 

It  a spaee  ot  piecewise  iv>l  ynonnu  1 s of  one  hiohoi  deuree  is  used  tor  the  pressure 

than  tot  the  concentrat  ion , optimal  ordet  1*  estimates  are  obtained.  Corollary  4.4 

- 1 .' 

■ tates  that  by  u(\latino  the  preeondi  t ioninq  matrix  l.  from  (A.1')  every  (At' 
time  stejs,  two  iterations  ot  the  equations  for  concent  rat  ion  will  suttioe  tv,'  obtain 
■ -ut  it  suit-;.  finally  extensions  tv'  methods  with  time  truncation  errors  ot  ordet 
0((At)fc)  'it  e vf  i soussed  brietlv. 

Theorem  4.1:  1 et  (c,p)  satisfy  (1.1)  — (l.S>  and  (Cl',l,n)  satisfy  (.1.20), 

(!..’»'),  and  t t.A“)  anil  let  k be  fixevi  as  in  C'..’9).  If  wo  obtain  norm  rediiv't  iv'ns 
in  ( t . .'4 ' of  t he  f orm 
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11  11  * ■ * i 
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where  D is  qiven  in  (4.36). 


1 roof : l.et  *.  = C - c and  n = P - p . We  shall  first  obtain  an  estimate 
on  Vn  at  the  At  -time  levels  in  terms  of  C and  the  error  induced  at  each  t iau> 
level  trom  t.he  iterative  approximation  to  (3.4).  Subtract  (2.10)  from  (3.20.b)  to 
obtain  at  each  At^-time  level, 
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f a . Since  C satisfies  (3.21)  . 


and  thus  (2.41,  we  can  use  (2.5), 


(2.14),  (2.15),  (3.22),  and  (3.23)  to  obtain 
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We  then  see  that 
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Wo  next  consider  the  concentrat ion  equations.  Recall  that  d r ^ 

t At 
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A--  in  in  older  to  most  efficiently  make  use  of  the  start.inq  proceduix's  i3.7.a 

of  the  con  lunate  qradient  iteration,  we  need  to  use  the  test  function 
\ ■ C Atd  c in  (4.5).  Then  we  see  that  with  a sliqht  rearranaomont , 

tin'  let t side  of  (4.5)  becomes 
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wheie  we  have  list'd  the  norms  and  semi-norms  defined  in  (3.23).  UsiiKi  techniques  like 
those  usexi  in  obtaininq  (4.41,  we  sum  the  first  two  terms  on  the  rioht  of  (4.5) 
from  n 0 to  n = C - l and  use  (2.14)  to  see  that 
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Next,  the  third  term  on  the  riqht  of  (4.5'  can  be  summed  and  split  as  follows: 
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X-4. 


We  note  that  from  (2.5.f)  and  (2.16),  we  can  bound  the  second  term  on  the  right  of 


(4.8)  by 
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We  shall  now  state  an  induction  hypothesis.  Assume  that  for  m - 0,1,..., If  - 1/k]  , 

(II) 


|ypmkll  . i 2k5 

L 


where  Ks  is  given  by  (2.15).  Clearly  from  (2.15)  and  (3.27)  we  know  that  (II)  is 


satisfied  for  m = 0.  Then  from  (2.28),  we  see  that  under  (II), 
(4.10)  l|£’uvplnk|l  < 4K  , m = 1,2,...,  [£-l/k 5 


for  use  in  (4.9).  We  next  bound  the  first  term  on  the  right  of  (4.8)  as  follows: 
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- utc.fijW*)  + u(cn  ,£^VP  ) - u(C  ^kVP  ) ) • Vc  ,dtc  )|At  = T3  + T4 

We  then  use  (2.14)  and  (2.16)  to  bound  the  second  term  on  the  right  of  (4.11).  Note 
that 
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We  next  use  summation  by  parts  to  treat  the  first  term  on  the  right  of  (4.11).  Note 
that 
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have 
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Wr  next  use  (2.3),  (2.12),  and  (2.14)  to  bound  the  first  term  on  the  right  of  (4.14), 
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The  third  and  fifth  terms  on  the  right  of  (4.14)  can  be  bounded  in  a similar  fashion 
to  yield 


(a(c1)V(p  - p) 1 • V (c  - c)1,^)  | + |(a(c%V(p  - pJ^-Vfc  - 
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k23{  II  !l  1 + h2r+2s“6}  . 


We  integrate  by  parts  and  use  (2.12)  and  (2.15)  to  bound  the  second  term  on  the  right 
of  (4.14).  Denote  a(cn)  by  a11  in  what  follows. 
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Note  that  by  (2.9) , Y a < K . In  order  to  bound  the  second  term  on  the  right  of 

l»n  - 5 


n=0 


(4.17)  we  shall  need  to  introduce  an  auxiliary  problem  to  perform  a Nitsche  lift  [14,10] 
2 

Let  i);  £ H (ft)  be  the  solution  of  the  elliptic  problem 
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(4.19) 


where 


(4.20) 


, . , n dip  n 

b)  a(c  } T7  = y ' 


x e Dft  , 


n n n-l__  n-1 

Y = a £ Vdfcc  "V 

= G1cn_1  e H1/2(3ft),  n = 1 , . . . ,N 


For  ,»  definition  of 


set*  [13],  Note  that  since 


1/2 


11  * (ihl)  and  its  norm  | • | , 


9c  2 + f 

y ‘ L (Jill  ( i- ) ) , we  have 

(4.21) 


°1 " *»  1/2+e  — ^25 

L (J,H  7 OSD) 


Then  under  the  smoothness  assumptions  on  a and  c,  we  use  (S)  and  Lemma  2.2 
to  obtain 
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w«'  then  use  V = (p  - p)r  * in  (4.23)  and  (2.10)  to  see  that  from  (4.17), 
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Similar  estimates  yield  the  bounds 
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Next,  for  the  vector-valued  function  u = u(x,c,q)  : il  x [0,1]  x R ->  ]R  , we  define 

2 1 2u.  . n 

a)  u*  (x)  l / ~r“-  (x,0cn  + (1  - 0)cn  , f^-)d0  , 

Un  i=i  o nc  3xi 


(4.27) 


2 1 2u.  _ -n 

u’  (X)  } / --*•  ( x , 0 c n + (1  - 0)cn_1,  V^—  ) d 0 , 

Z,\\  • i „ l’C  (’X. 


•)  — a(o),  i = 1,2  . 

3qi 


Usinq  (4.27),  we  treat  the  second  term  in  (4.13)  as  follows: 
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Then  using  the  same  techniques  as  above,  if 
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We  next  obtain  bounds  for  the  last  two  terms  on  the  right  side  of  (4.5)  using  (3.24) 
Since  different  starting  procedures  were  used  in  the  conjugate  gradient  iteration 
to  obtain  and  C°  for  n ^ 2,  we  shall  estimate  each  case  separately.  From 
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2 ii  r 2 nii  , , ii  , .nn  . ....2m  , .m 


(4.34) 


1 p;t||«  Cn||  ♦ (Atr  ||«V1||  nH  ||dtcn||v?  + (At)^||dt;n||  n} 

b b 

1 ^ j 
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+ At 
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* <At)2||dt;f.n|| 2n)  + (at) 
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at 


oo  1 
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Combining  the  above  estimates,  we  obtain 

J T Hv"ii>  * i *jl  ‘IIs 


i -»»v  . i’j‘  n|C"*iii2„-  iit“n2 ) 

n=0  ~ * n=0  b b 


< 1,  ii  _ 4 n 2 

-8  b 


4-1 


+ k4o  i <-lfB^t>ii'nii; 

w n=l 


(4.35) 


+ K4!{  II  C°  II  f + h2r  + h2S  + h2r+2s_(’  + (At)2  + (At1)4D*2) 


H-1AD 


♦ " i " iiv<pmk  - »™k>r'At  , 
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m=0 


where  in  this  case 

(4.36) 


D 


*2 


x /llcll 


2 1 ' " ~ " 2 1 
W ( J; H ) W (J;H  ) 


is  to  be  used  in  (2.29)  to  determine  k qiven  At.  We  next  consider  the  precondi- 
tioned conjuqate  gradient  iteration  on  the  pressure  equations  to  treat  the  last  term 
on  the  right  of  (4.35).  From  (3.23),  (3.24),  and  (4.4),  with  m ■*  2,  we  see  that 


-28- 


||v(?"k-P"kl|l  i a.2  Ilf"*-  r*kl|  ^ i*.2e2l|f”'k-  2r("-1,k.r'”-2lk| 
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£ 

* 2 f n mk 

< (a  /at)  pA  H'l 


I „k  ' -'ll"lm'nk|l  ta-„K  * 

a a 


(m-2)  k I 


(m-2) k 


+ (At^' 


1J1 

3t2 


o»  X 
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} * a* 


2P2l|F"k-  r"k||  mk 

a 


(4.37) 


iKal(||t"k||  * lie 


4 J 


(m-Uk||  + ||cC-2)k|(  + (Ati>2  + hr} 


1 

♦ 2lr^-]2  Ml  r'mK  - pmKH 


mk 


+ 4 


t *\2 

|a_  II  T <m~ 1 ) k p(m_1)k|| 

a 11  11  (m-l)k 


t 


♦ 2 


” (m-2) k _ ^ (m-2) k i 


(m-2) k 


Wo  thou  set  that  if 


(4.38) 

then 


1 *2 

P2  - 12”  <a*/a  ) > m ^ 2 , 


af-l/kl) 


(4.39) 


l/KJI  ml,  2 f i-1  2 ~ i 

l Atll?n^||  <k  { l At || cn||  + h2r  + (At  )4n‘l. 

= 2 **  ln=0  1 > 


A similar  result  will  hold  for  m = 1 if  we  iterate  sufficiently  many  times  that 


(4.40) 


P2  (At)  , 


m = 1 . 


.1 


In  order  to  apply  a discrete  version  of  the  Gronwall  Lemma  in  the  H norm,  we  must 
first  shift  the  coefficients  in  the  third  term  on  the  right  of  (4.e)  to  obtain  a 


telescoping  sum  in  the  ||  • ||  semi-norm  and  then  introduce  an  L“  term  into  the 
telescopinq  sum  to  obtain  H terms.  Note  that 


-2a- 


(4.41) 


|cn||2  ,+  ((b(tn+1)  - b(tn)jvcn,vtn) 

1 ",  n-1 


i !Unll2n_i  + WWfK 

b 


I'll  2 


Also  note  that 


(4.42) 


iu"'1||2  - III "II2  - 2At(dtt.".t")  . iat)2||dtt;"ir’ 

i&ll-l'-ltf  * A«16(Un||2  • 

We  now  sum  (4.41)  and  (4.42)  from  n = 0 to  n ■ t,  add  the  results  to  (4.35),  use 
(3.26),  (4.3°)  and  (4.40),  and  apply  a version  of  the  discrete  Gronwall  Lemma  in  the 
|| . ||  ♦ || . ||  norm  (equivalent  to  the  H norm)  , to  obtain 


(4.43) 


<-l 

l 

n=0 


l (At ||d  tn || 2 + (At)2||d  Cn||2  ) + lie'll^ 

_ t.  ^ b*1 


c II  2 


2r 


2s  . 2r*2s-6 


. 4*2 , 


KAnlh"  + h‘“  + h ♦ (At)  * •'  (At  ) D ) . 
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Since  (4.43)  holds  for  any  « = 0,...,N,  wo  can  replace  II C. ? II T by  sup  ||  C ||.- 
Then  we  see  from  (2.3),  (2.15),  (4.4),  (4.3/)  and  (4.43)  that 


0<n<N 


|v,'»k||  II  Vnmk  l|  IhPII  . „ IS* 

L L b (0 : b ) 


(4.44) 


< K . K„h->K(hr  * lAt.iV  * |k”kll  • IK(”,’l'kll  * Ik 

— SO  I 


(m-2) k , 


< K + Kh_1(hr  + hS  t At  t (At  )2D*1 
— S l 


Then  if  r,s  > 2,  At  = (At  >"'P*,  and  At  • h' for  f ■*  0,  we  see  that  for  At 

and  h sufficiently  small  our  induction  hypothesis  will  be  satisfied  for  all 

m < jN/kJ.  Finally  (2.14),  (4.43),  and  the  tr innate  inequality  yield  (4.1),  the  desired 

result . 

We  note  that  if  we  replaced  in  the  coefficient  u in  (3. 20. a)  by 

vp|n/kfl  or  VF.l'®1’  k®  (as  in  (.'..'/)  and  the  discussion  followinq  (2.27))  wc  would 

• i , * 2 2 

replace  the  trim  It  (At  1 in  (4.1)  by  P (At ^ 1 whcic  P 


is  si  iqht  l y >1  i I t oi  out 


than  D . However,  the  use  of  the  preconditioned  conjugate  gradient  itei  at.  ion  with 


this  modification  would  still  involve  bounds  on  ^ 

3t ' 


HI 


>»  1 
L ( J ; H 


I.  ( J ; H ) 


Since  these  terms  would  probably  dominate  both  D and  P , t lie  ability  to  choose 
a much  larger  k and  thus  solve  (3.20.b)  much  less  frequently  would  motivate  the 
use  of  the  different  extrapolations  defined  in  (2.28).  However,  tin'  repeated  evalua- 
tion of  the  coefficients  u (and  b in  Case  11)  is  often  computational  1 y expensive . 

By  considering  the  u appearing  in  (l.l)  and  (1.2)  as  a separate  variable  to  be 
determined  from  the  elliptic  equation  (1.2)  by  a mixed  method,  we  could  I'ossibly 
extrapolate  the  variable  u instead  ot  just  Vp  to  achieve  both  tin'  accuracy  and 
the  computational  efficiency  of  not  having  to  evaluate  a new  coefficient  at  each  At- 
tune step.  This  idea  is  being  considered  elsewhere. 

We  note  that  if  b = b(x)  only,  then  the  l.n  defined  in  (3.3)  will  be  independent 
of  t and  n.  In  this  case  1.  i t AtR  will  be  factored  and  the  preconditioned 
conjugate  gradient  iteration  for  the  parabolic  equation  will  be  replaced  at  each  step 
by  a simple  back-solve  of  the  factored  1..  The  iteration  will  still  be  used  for  t tie 
elliptic  equations.  We  obtain  the  following  simplification  of  Theorem  4.1. 

Corollary  4.2:  Assume  b b(x)  and  k is  fixed  as  in  (2.2'').  It  we  obtain 
norm  reductions  in  (3.24)  (c)  and  (d)  of  the  form 


a)  p_  v (,\t  ) for  )'  , and 

" 1 


1 '*  r em_k 

h)  - 77  v lor  1 


then  there  exist  positive  constants  K,„  = K.  , (\,a.  ,b.  . ,K.  , i s 10),  h . and  i , 

4R  48  * * * t 0 0 

such  that,  if  At  *•  t and  h «•  h , 

0 0 


sup  < ||C  - c||*  • h II C - c||;  ) I 


(C  - .V'||\M 


K , i l,‘  * ' li‘S  * 1 ' (At)-  > n*‘  (At  t ) 4 ) 


(4 .4M 


who  r e 


2 l /Hell  i , 

W (J;H  ) W (J;H  ) 


We  next  consider  the  version  of  (4.5)  for  Case  II,  where  b = b(x,t)  is  replac'd 
by  b » b(x,c,Vp).  in  tins  case  (4.5)  is  replaced  by 


(4.4e) 


r -i  n f 1 

«dtin.x>  * (bic  ",ej’vi’")vi"*',vxi  - M Vjej  * 
• (|b(cn*,.V|> - b(C*",£”vp")]Vc‘*1.Vx) 

. - urc'VJVVvr-.xi 

. W . xl  * (b(C*n,^7Pn)7(cntl  - V»l, 


,n+l 

•) 


* aF — * x>  f (b(c  '\VF  >*<c  - C ).VX).  V t \ • 

A test  function  of  the  form  r.1'  ^ - £n  is  required  to  make  most  efficient  use  of  the 
iterative  procedure  defined  in  Section  2.  However,  this  choice  of  test  function 
causes  serious  problems  in  the  treatment  of  the  third  term  on  t he  riqht  of  < 4 . 4 «.■» ) . 
Standard  techniques  for  treatinq  a term  of  this  type  like  summation  by  parts  in  time 
(see  [b,9])  will  not  work  due  to  the  inability  to  treat  the  resulting  terms  of  the  foim 


r . , n+1  n-t-l  n+1  n,,v_n.  . n n.  . n __v  n-i  . -n  n 

l (lb(c  ,Vp  ) - b (c  ,VP  P ) - (b  (c  ,Vp  ) - b (c  .V.^P  ) ) ) Vc  ,VdtC,  ) At 

n=0 

For  this  reason,  in  the  proof  of  the  following  result,  we  shall  use  a test  t unction 
of  the  form  x = £n+1. 

Theorem  4.3:  Let  b = b(x,c,Vp)  and  (C^.P0)  satisfy  (3.20)  (with  b*'  1 


replaced  by  b(C  .f^VF  )),  (3.2b)  (with  |j  • ||  ^ replaced  by 


, and  (3.27)  and 


let  k be  fixed  as  in  (2.29).  If  we  obtain  norm  reductions  in  (3.24)  of  the  form 


a)  p'  < 7 At  for  n " 1 , 

1-4  - 


b)  p < (At) “ for  P , 
2 1 


I f-*l 


c)  p?  v r?  r*  for  p 


then  there  exist  iwsitive  constants  K.„  = K (X.a  , b ,<f  ,K.,i  < 10)  such  that,  if 

49  49  1 


At  ^ and  h <_  h^. 


(4 


.47)  sup  (||c-  c||2  + h||c-  c||2}  I K4i>th“,r  + h‘S_2  + (At)2  + (At^D*2}  , 


where  D is  given  by  (4.36). 

Proof : We  shall  use  the  same  notation  as  in  the  proof  of  Theorem  4.1.  We  obtain 

(4.4)  just  as  before  and  then  consider  (4.4b).  With  X = An+* , the  left  hand  side 
of  (4.46)  becomes 

t (IU"*lIL  - llc”IU  * IU"+ll|2„ 


2 At 


^4.481 


= (V'dtcn,i;,1+1)  + (b(c*n,^vrn)Vcn+1,vc.n+1)  . 
We  bound  the  first  and  second  terms  on  the  right  of  (4.46)  as  follows 


(4.4a) 


r 

r ^ n+1  1 

3c  , ' n 

, - , n+1  n+l) 

* 

L at  - V j 

+ \ (c  - c)  , C 

i s«{*2r 


3 c 

at2 


(At) 
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We  then  bound  the  third  and  fourth  terms  on  the  right  of  <4. 4b)  by 

I . . n + l . n + 1 *n  „v  n , „~n+l  „ n+1  n+1  m l -n+1 

| ( lb  (c  ,Vp  ) - b (C  ,B.VP  )]Vc  ,VC  ) + ([u(c  , Vp  ) • W 

, *n  ^v.-^n.  „ n,  n+l,i 
- u(C  »Ekvr  I’Ve  l,r,  ) j 


(4.50) 
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+ ( At.  j ) 


Lj 
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+ irl  + h"s  ‘ + ||t.nl|  + 1 • 


ao  1 

L (J ; H ) 


Note  we  have  used  our  induction  hypothesis  (II)  in  (4.50)  as  before.  We  next  use 
to  bound  the  last  two  terms  on  the  right  of  (4.4b).  For  n ' 0,  we  multiply  by  it 
and  obtain 


- X x- 
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- ( ~C  , n+ll  , ,*n  .A'  ti...  ,n*l  -n+1,  „ n + 1 
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At 
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Then  »t  wo  Herat  c suf  f ioiont  l y many  times  that 
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(4.S.') 

wo  soo  that  is  bounded  by 
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ln7  ilknM 
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‘Atf  H^ll'nMb-  * V''lk"!L  * . 
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Wo  then  sum  the  above  inequalities  on  n,  oombino  the  results,  and  lot  P be  as 
in  (4. 4t>)  to  obt  at  n 

(i->..lk‘ii;.  T„|y -{u.^}Jt)|ir1ir'„ 

n~  0 v v * b 

r .... , 4,  , c-i 
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The  rest  of  the  proof  follows  as  i ti  the  proof  of  Theorem  4.1. 


1 


We  note  that  tsee  |t-|'  it  we  updat  e the  pi  eeond  1 1 tom  no  matrix,  oaeh  (At)  * -t  inu 
steps,  we  obtain  a norm  reduet  ion  of  P(»At)  with  one  iteration  and  (''tAtl  with  two. 
Usitw  this  idea  we  obtain  the  tollowmq  eoiollaiv  to  Thtn'i  em  4.'. 

Oorol  I arv  4.4:  be*  tv'  ,P  ' -atistv  (i..'A),  t 5 . .'el  and  t*..'-'  as  in  Vheoi  em  4.'. 


by  updat  i u a ! 


1*  at  oaoh  i it 


1 


Vhool  em  l.t  with  onlv  two  ifoiat 


•time  tops,  we  obtain  hyjxM  lies  i s a'  ot 
,<.'*)  el  t line  st  el- . 


1 1 


Using  techniques  like  those  above  wo  can  obtain  corresjoiiii  mg  lesults  for  the 


Crank-Nicolson-Galerkin  approximations  to  (1.1)  — (1.5)  . These  methods  do  not  have  the 
;aine  stability  properties  as  the  methods  analyzed  above  unless  we  make  a restriction  like 

(4.55)  At  c k h2  . 

~ 54 

Under  the  assumption  (4.55)  and  slightly  stronger  smoothness  assumptions  on  p and  c, 

2 4 

we  would  obtain  the  same  types  of  results  as  above  with  (At)“  replaced  by  (At)  . 
Without  (4.55),  we  would  need  to  iterate  more  to  obtain  the  necessary  norm  reductions 
from  (3.24).  See  (6]  for  particulars. 


5.  Computation.!  1 Considerations 

In  this  section  we  shall  consider  some  rough  operation  counts  to  estimate  the 
computational  complexity  of  the  methods  presented  here.  We  shall  see  that  tin  ii  Midi 
tinned  conjugate  jradient  iterative  methods  presented  allow  us  to  obtain  neat  ojitini.il 
order  work,  estimates.  Thus  these  methods  are  very  efficient  computationally. 

Recall  that  we  have  two  space  variables  (d  - 2).  George  (111  has  shown  in  a*  i.ie 
special  cases  that,  with  M = M (hi  - dim  A*  , the  procedure  of  setting  up  and  factoring 

3/2 

4 + * t.B^  requires  0(M^  ) operations  and  that  the  solution  of  (3.3),  givi-n  tin. 

factorization,  requires  0 (M ^ log  M^)  operations.  Similarly,  the  work  involved  in 

setting  up  and  factoring  A and  solving  (3.4!  using  this  factorization  are  0(M  ") 

n 2 

and  OIK,  log  M ( ) respectively.  Hoffman,  Martin  and  Rose  112]  have  shown  that  such 

bounds  are  minimal.  Thus,  if  we  conjecture  the  validity  of  the  above  estimates  for 

our  jroblt'm  and  refactor  4 t AtB  and  A and  solve  (3.3)  and  (3.4)  at  each  time 

n n 

step , the  total  amount  of  work  done  using  the  method  presented  in  [10]  is 

1 2 3/2  3/2  3/2 

(5.1)  0(Nt((M1  + Mj  log  M ^ ) + (M.?  + M.,  loa  M.,)})  = OIN^.Jm^  + M./  ))  , 

1 

wnere  NT  is  the  total  number  of  time  steps  (N^,  * ) . We  note  that  the  work  of 

refactor izat ion  dominates  the  estimates. 

Using  the  preconditioned  conjugate  gradient  iterative  method  presented  here,  one 
does  not  need  to  refactor  at  every  time  step.  Instead  only  one  factorization  of 


L • + AtB  and  A need  be  done.  Also,  using  the  different  time  increments  for 
VJ  U 0 

pressure  and  concentration,  we  need  only  solve  (3.4)  at  every  k*'*1  time  step;  thus 
the  total  number  of  times  (3.4)  need  be  solved  is  N^yk . Thus  letting  tc  and  k 
be  the  number  of  pro-condi tioned  conjugate  gradient  iterations  needed  to  achieve  the 


norm  reductions  given  in  (4.0)  (k  and  t n are  constants,  independent  of  h,  n 
and  At) , we  see  that  the  total  work  required  for  Case  I is 


Then  since  ^ ^ 


1 r/2 

— = 0(Mj  ),  we  see  that  even  for  r = 2 (piecewise  linear 


elements),  the  second  and  fourth  terms  in  (5.2)  dominate  and  the  work  is 

N_ 


(5.3) 


°(VlMi  loq  M]  + — k'2M0  loq  M2) 


We  note  that  the  number  of  unknowns  in  the  problem  is 

N 

(^••1)  o(ntmi  + y-  mj  , 

so  (5.3)  represents  nearly  best  possible  order  work  estimates. 

It  is  computationally  wasteful  to  iterate  exactly  times  at  each  time  step 

(respectively  times  at  every  k time  step)  in  order  to  achieve  the  pessimistic 

bounds  qiven  in  (4.0).  Instead,  one  can  monitor  the  norm  reduction  actually  produced 
at  each  step  of  the  iteration  and  stop  iterating  when  sufficient  norm  reduction  is 
achieved.  Additional  stopping  criteria  can  be  imposed  in  this  monitoring  process. 


See  [f>)  for  a discussion  stopping  criteria  for  related  problems. 
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